function [Rspec_sig, freq, Rspec, Ispec,v2] = texture_spec(fname, dflag, mmax)
%[Rspec_sig, freq, Rspec, Ispec,v2]=texture_spec(fname,dflag, mmax) deadleaves NPS
% Computes Dead Leaves noise-power spectrum from image file for deal leaves
% test chart WITH NOISE CORRECTION. The following operations are included,
% ROI selection, 2D NPS estimation, radial NPS integration.
%
% fname = path for image file name OR array of (nxn) image pixel data.
%         If a color image is chosen, a luminance color report will be
%         computed and used.
% dflag = 0 no 2D detrending of data array before NPS estimation
%       = 1 detrend
% mnax  = (optional)size of data array used (mmax, mmax) used for NPS
%         estimation, default = 256
% Rspec_sig = 1D (radial) signal spectrum, corrected for noise spectrum
% freq  = spatial frequencies (vector) for sampling of Rspec_sig
% Rspec =  1D (redial) signal spectrum, NOT corrected for noise spectrum
% Ispec =  2D uncorrected noise-power spectrum
% v2 = variance computed fron NPS, (used for debugging spectrum estimates)
%
%Needs: imageread, getroi, detrend2, deadleavesNPS
%
% Peter D. Burns, pdburns@ieee.org 5 Nov. 2015

if nargin<3;
    mmax = 256;
end

if nargin<2
    dflag = 0;
end
if ischar(fname)==1
  [status, dat] = imageread(fname);
  dat2 = getroi(dat,'Select deadleaves texture region');
end
dat2 = double(dat2);

[nlin, npix, nc] = size(dat2);
if mmax~=0
    % Simple optional data cropping
    nlin = min(nlin,mmax);
    npix = min(npix,mmax);
    dat2 = double(dat2(1:npix,1:nlin));
end
if nc == 3;
    dat2 = rgb2lum(dat2);
end

if dflag==1
  % Simple 2D linear (a plane) subtraction
  dat2 = detrend2(dat2, 1, 1, 0);
  disp('detrending signal')
end
[Rspec, freq, Ispec, v2] = deadleavesNPS(dat2, 1, 0);

dat3 = getroi(dat, 'Select uniform noise region');
dat3 = getroi(dat3,'Refine noise selection (if needed)');
dat3 = double(dat3);

if nc == 3;
    dat3 = rgb2lum(dat3);
end

if dflag==1
  dat3 = detrend2(dat3, 1, 1, 0);
  disp('detrending noise')
end
[Rspec1, freq1] = deadleavesNPS(dat3, 1, 0);

% Interpolate noise spectrum to same size as signal spectrum
Rspec1int = interp1(freq1, Rspec1, freq, 'spline');
Rspec1int = clip(Rspec1int,eps,inf);

% Subtract noise spectrum
Rspec_sig = Rspec - Rspec1int;
% Avoid negative spectrum values
Rspec_sig = clip(Rspec_sig,eps,inf);

%Plotting of results
% figure,
% semilogy(freq, Rspec1int,'r--','LineWidth',1.6), hold on
% semilogy(freq, Rspec,'LineWidth',1.6)
% semilogy(freq(1:2:end), Rspec_sig(1:2:end),'k:','LineWidth',1.6)    
% hold off
% mmax = 1.05*max(Rspec);
% mmin = mmax*1e-6;
% axis([0 .5 mmin mmax])
% 
% xlabel('Frequency, cy/pixel')
% ylabel('Power Spectrum')
% legend('Noise', 'Texture signal','Corrected')
